Lines & Polygons

 

Rodney Dyer, PhD

Learning Objectives

This section covers the following topics:

  • Creation of sf::LINESTRING and sf::MULTIPOLYGON objects from raw data.

  • Loading map polygons from get_data()

  • Importing data from ESRI shapefile.

  • Implementing Spatial Joins.

Lines

Data Sources: De novo Creation

It is possible to make data from scratch by entering coordinates directly.

df <- data.frame( id = c(rep("Rodney",5), rep("Sarah",5)),
                  Data = rnorm( 10, 42, 42),
                  Longitude = rnorm(10, -78, 1 ),
                  Latitude = rnorm(10, 37, 1) )
df
       id        Data Longitude Latitude
1  Rodney -46.1064984 -78.77667 37.50880
2  Rodney  56.4871453 -79.14573 36.87373
3  Rodney  14.0225502 -76.88675 39.18777
4  Rodney 138.5959504 -78.58731 37.84106
5  Rodney  27.0122645 -77.52085 37.67201
6   Sarah   0.9570969 -78.21445 35.69324
7   Sarah   5.2697847 -78.61477 37.39748
8   Sarah  70.0049514 -78.26441 36.21965
9   Sarah   7.5430525 -79.16548 37.08932
10  Sarah  42.7788551 -77.18989 36.76554

Data Sources: De novo Creation

To Make to sf as POINT objects as normal.

df %>%
  st_as_sf( coords = c("Longitude", "Latitude"), crs=4326 ) -> pts
pts
Simple feature collection with 10 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -79.16548 ymin: 35.69324 xmax: -76.88675 ymax: 39.18777
Geodetic CRS:  WGS 84
       id        Data                   geometry
1  Rodney -46.1064984  POINT (-78.77667 37.5088)
2  Rodney  56.4871453 POINT (-79.14573 36.87373)
3  Rodney  14.0225502 POINT (-76.88675 39.18777)
4  Rodney 138.5959504 POINT (-78.58731 37.84106)
5  Rodney  27.0122645 POINT (-77.52085 37.67201)
6   Sarah   0.9570969 POINT (-78.21445 35.69324)
7   Sarah   5.2697847 POINT (-78.61477 37.39748)
8   Sarah  70.0049514 POINT (-78.26441 36.21965)
9   Sarah   7.5430525 POINT (-79.16548 37.08932)
10  Sarah  42.7788551 POINT (-77.18989 36.76554)

Plotting Points

We can easily mix individual aesthetic properties to change the representaion of the points provided by using geom_sf().

 

ggplot( pts ) + 
  geom_sf( aes(color = id, size=Data)) + 
  coord_sf() + 
  theme(axis.text.x = element_text(angle = 45, 
                                   vjust = 0.5, 
                                   hjust=1))

 

 

Data Sources: De novo Creation

To make sf as LINESTRING we can group by the id value and then we need to summarize the Data component, then cast it to a Linestring.

df %>%
  st_as_sf( coords = c("Longitude", "Latitude"), crs=4326 ) %>%
  group_by( id ) %>%
  summarize( Data = mean( Data ) ) %>%
  st_cast("LINESTRING") -> lines

 

lines
Simple feature collection with 2 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -79.16548 ymin: 35.69324 xmax: -76.88675 ymax: 39.18777
Geodetic CRS:  WGS 84
# A tibble: 2 × 3
  id      Data                                                          geometry
  <chr>  <dbl>                                                  <LINESTRING [°]>
1 Rodney  38.0 (-77.52085 37.67201, -79.14573 36.87373, -78.77667 37.5088, -78.…
2 Sarah   25.3 (-78.21445 35.69324, -78.26441 36.21965, -77.18989 36.76554, -78…

Visualizing Line Objects

Just like when we were working with POINT objects, we can use both built-in graphing approaches as well as ggplot() for visualizing. Here I’m goint to plot the two line objects using normal graphics.

plot( lines["Data"], 
      col=c("red","green"),
      lwd = 2)

 

 

Operations on LINESTRING Objects

Physical length of the LINESTRING objects

st_length( lines )
Units: [m]
[1] 497817.9 376535.9

The bounding box around all the lines.

st_bbox( lines )
     xmin      ymin      xmax      ymax 
-79.16548  35.69324 -76.88675  39.18777

The area of the convex hull encompassing all the lines.

library( units )
lines %>% 
  st_union() %>% 
  st_convex_hull() %>% 
  st_area() %>%
  set_units( km^2 )
39386.81 [km^2]

Textual Versions of geometry

While it is easy to go from text \(\to\) POLYGON, the same thing can be done going in the opposite direction.

lines$geometry[1] %>% st_as_text()
[1] "LINESTRING (-77.52085 37.67201, -79.14573 36.87373, -78.77667 37.5088, -78.58731 37.84106, -76.88675 39.18777)"

Polygons

Polygons

Polygons are simply lines whose first and last coordinate are exactly the same

df.p <- df[ c(1:5,1, 6:10,6),]
df.p
        id        Data Longitude Latitude
1   Rodney -46.1064984 -78.77667 37.50880
2   Rodney  56.4871453 -79.14573 36.87373
3   Rodney  14.0225502 -76.88675 39.18777
4   Rodney 138.5959504 -78.58731 37.84106
5   Rodney  27.0122645 -77.52085 37.67201
1.1 Rodney -46.1064984 -78.77667 37.50880
6    Sarah   0.9570969 -78.21445 35.69324
7    Sarah   5.2697847 -78.61477 37.39748
8    Sarah  70.0049514 -78.26441 36.21965
9    Sarah   7.5430525 -79.16548 37.08932
10   Sarah  42.7788551 -77.18989 36.76554
6.1  Sarah   0.9570969 -78.21445 35.69324

Polygon de novo Creation

df %>%
  st_as_sf( coords = c("Longitude", "Latitude"), crs=4326 ) %>%
  group_by( id ) %>%
  summarize( Data = mean( Data ) ) %>%
  st_cast("POLYGON") -> polygons

 

polygons
Simple feature collection with 2 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -79.16548 ymin: 35.69324 xmax: -76.88675 ymax: 39.18777
Geodetic CRS:  WGS 84
# A tibble: 2 × 3
  id      Data                                                          geometry
  <chr>  <dbl>                                                     <POLYGON [°]>
1 Rodney  38.0 ((-77.52085 37.67201, -79.14573 36.87373, -78.77667 37.5088, -78…
2 Sarah   25.3 ((-78.21445 35.69324, -78.26441 36.21965, -77.18989 36.76554, -7…

Visualizing Polygon Objects

ggplot( polygons ) + 
  geom_sf( aes( fill = Data ), alpha=0.75 ) + coord_sf() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1))

Polygons From Data Frames

map_data( "county", "virginia") %>%
  dplyr::select( Longitude = long,
                 Latitude = lat,
                 group,
                 County = subregion) -> va_counties
head( va_counties )
  Longitude Latitude group   County
1 -75.27519 38.03867     1 accomack
2 -75.21790 38.02721     1 accomack
3 -75.21790 38.02721     1 accomack
4 -75.24655 37.99283     1 accomack
5 -75.30384 37.94127     1 accomack
6 -75.31530 37.92981     1 accomack

ggplot( va_counties, aes( Longitude, Latitude) ) + 
  geom_point( size=0.25 ) + 
  coord_quickmap()

va_counties %>%
  filter( County %in%  c("hanover","henrico") ) %>%
  ggplot( aes(Longitude, Latitude) ) + 
  geom_polygon( aes( fill = County), alpha=0.1 ) +
  geom_point( aes( color = County) ) +
  coord_quickmap()

ggplot( va_counties, aes( Longitude, Latitude) ) + 
  geom_polygon( aes( group=group ), fill="grey80",
                color = "black", size = 0.25) + 
  coord_quickmap()

ggplot() + 
  geom_polygon( aes( Longitude, Latitude, group=group ),
                fill="grey80",color = "black", 
                size = 0.25, data=va_counties) +
  geom_sf( aes( color=Data), 
           lwd=3, data=lines) + 
  coord_sf()

Shapefiles 🤷

 

Shapefiles

The ESRI Shapefile is a non-standardized format for packaging vector based data (POINT, LINESTRING, POLYGON, etc.).

However, it is not actually a file but it is a collection of files, which may (or may not) expand directly in the folder or within a subfolder.

Archived ZIP Files.

I’ve uploaded some shapefiles to the Github Repository for this class. These are from the Richmond City GIS Department and represent the centerlines of roads in the city as well as the zoning districts.

Here are the URL’s for both of these files.

roads_url <- "https://github.com/dyerlab/ENVS-Lectures/raw/master/data/Centerlines-shp.zip"
district_url <- "https://github.com/dyerlab/ENVS-Lectures/raw/master/data/Zoning_Districts-shp.zip"

Downloading Files Using R

You can use your browser and copy and paste those addresses in to download the files to your computer. Then you’ll have to move those archives to your PROJECT FOLDER (you are using Projects right?).

OR

You can have R download the file and save it locally.

download.file( district_url , destfile = "./Districts.zip")
download.file( roads_url, destfile =  "./Roads.zip")

Unziping Files

You can open your computer finder and have the OS unzip the archives.

OR

You can have R do it.

unzip("Districts.zip")
unzip("Roads.zip")

File Components

Loading A Shapefile

You can load it in using sf::st_read() and passing it the .shp file path to it.

st_read("Zoning_Districts.shp") -> districts 
Reading layer `Zoning_Districts' from data source 
  `/Users/rodney/Public/github/teaching/Shapefiles/Zoning_Districts.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 634 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 11743500 ymin: 3687944 xmax: 11806060 ymax: 3744741
Projected CRS: NAD83 / Virginia South (ftUS)

Shapefiles Convert to sf Objects

By default, the data associated with a vector object are put into a data.frame and the geometry object is properly created.

head( districts, n=2 )
Simple feature collection with 2 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 11773600 ymin: 3730159 xmax: 11789510 ymax: 3731016
Projected CRS: NAD83 / Virginia South (ftUS)
  OBJECTID Name Ordinance OrdinanceP Conditiona AdoptionDa Comment
1        1 RO-2      <NA>       <NA>         No 2000-01-01    <NA>
2        2  B-2      <NA>       <NA>         No 2000-01-01    <NA>
           CreatedBy CreatedDat             EditBy   EditDate
1 richard.morton_cor 2020-08-24 richard.morton_cor 2020-08-24
2 richard.morton_cor 2020-08-24 richard.morton_cor 2020-08-24
                              GlobalID Shape__Are Shape__Len
1 334799f0-fe38-46bf-97c2-260f5a036559   60150.29   983.6815
2 558df9cd-4f9c-4248-a689-bc2d9c79d060   56987.01   971.8832
                        geometry
1 MULTIPOLYGON (((11773598 37...
2 MULTIPOLYGON (((11789222 37...

Simplifying Data

names( districts )
 [1] "OBJECTID"   "Name"       "Ordinance"  "OrdinanceP" "Conditiona"
 [6] "AdoptionDa" "Comment"    "CreatedBy"  "CreatedDat" "EditBy"    
[11] "EditDate"   "GlobalID"   "Shape__Are" "Shape__Len" "geometry"

Simplifying Data

districts %>% 
  dplyr::select( -Comment, -CreatedBy, -CreatedDat, -EditBy, -EditDate ) %>%
  dplyr::select( -Shape__Are, -Shape__Len) -> districts
head( districts, n=1)
Simple feature collection with 1 feature and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 11773600 ymin: 3730159 xmax: 11773940 ymax: 3730493
Projected CRS: NAD83 / Virginia South (ftUS)
  OBJECTID Name Ordinance OrdinanceP Conditiona AdoptionDa
1        1 RO-2      <NA>       <NA>         No 2000-01-01
                              GlobalID                       geometry
1 334799f0-fe38-46bf-97c2-260f5a036559 MULTIPOLYGON (((11773598 37...

District Codes

districts$Name %>% unique() %>% sort()
 [1] "B-1"     "B-2"     "B-3"     "B-4"     "B-5"     "B-6"     "B-7"    
 [8] "CM"      "DCC"     "I"       "M-1"     "M-2"     "OS"      "R-1"    
[15] "R-2"     "R-3"     "R-4"     "R-43"    "R-48"    "R-5"     "R-53"   
[22] "R-6"     "R-63"    "R-7"     "R-73"    "R-8"     "R-MH"    "RF-1"   
[29] "RF-2"    "RO-1"    "RO-2"    "RO-3"    "RO2-PE2" "RO2-PE4" "RP"     
[36] "TOD-1"   "UB"      "UB-2"    "UB-PE2"  "UB-PE3"  "UB-PE4"  "UB-PE5" 
[43] "UB-PE6"  "UB-PE7"  "UB-PE8"  "UB-PO1"  "UB-PO2"  "UB-PO3"  "UB-PO4" 
[50] "UB2-PE8"

Visualizing

plot( districts["Name"] )

Visualizing Overlays

plot( st_geometry( districts ))
plot( st_geometry( districts[districts$OBJECTID==530,] ), 
      col='red', add=TRUE )

Vector Operations

Spatial Joins

In the topic covering [Relational Operations], we used Primary and Foreign Keys to join data.frame objects and combine data. We can do simliar operations using geospatial positions to perform spatial joins.

There are a wide variety of operations available:

Check out the Cheetsheet

Secondary Vector Data Set

road.shapefile <- st_read("Centerlines-shp/tran_Carriageway.shp")
Reading layer `tran_Carriageway' from data source 
  `/Users/rodney/Public/github/teaching/Shapefiles/Centerlines-shp/tran_Carriageway.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 29081 features and 25 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 11734060 ymin: 3682790 xmax: 11817490 ymax: 3751927
Projected CRS: NAD83 / Virginia South (ftUS)
names( road.shapefile )
 [1] "FID"        "Carriagewa" "AssetID"    "StreetType" "Functional"
 [6] "FIPS"       "LeftFromAd" "LeftToAddr" "RightFromA" "RightToAdd"
[11] "PrefixDire" "ProperName" "SuffixType" "SuffixDire" "FullName"  
[16] "RouteName"  "OneWay"     "PostedSpee" "CADRouteSp" "CreatedBy" 
[21] "CreatedDat" "EditBy"     "EditDate"   "GlobalID"   "SHAPE_Leng"
[26] "geometry"

Check Projections

Dyer’s First Rule: Make Sure Projections Are the Same

st_crs( road.shapefile ) == st_crs( districts )
[1] TRUE

If they differed, one could do something like:

road.shapefile %>%
  st_transform( crs = st_crs( districts ) ) -> road.shapefile

Data Cleanup

road.shapefile %>%
  dplyr::select( FullName, OneWay, StreetType,
                 SpeedLimit = PostedSpee, Length = SHAPE_Leng,
                 geometry) %>%
  mutate( OneWay = factor( OneWay ),
          StreetType = factor( StreetType) ) -> roads
summary( roads )
   FullName          OneWay                    StreetType      SpeedLimit   
 Length:29081       1   :    2   Alley              : 3139   Min.   : 0.00  
 Class :character   FT  : 3225   Artery             : 5201   1st Qu.:25.00  
 Mode  :character   TF  : 3148   Highway            :  384   Median :25.00  
                    NA's:22706   Highway Interchange:   93   Mean   :25.17  
                                 Private            : 1540   3rd Qu.:25.00  
                                 Ramp               :  378   Max.   :55.00  
                                 Secondary          :18346   NA's   :7959   
     Length                   geometry    
 Min.   :    2.943   LINESTRING   :29081  
 1st Qu.:  164.849   epsg:2284    :    0  
 Median :  292.517   +proj=lcc ...:    0  
 Mean   :  379.338                        
 3rd Qu.:  461.146                        
 Max.   :17851.326

Visualizing: Filter Highways & Plot

roads %>%
  filter( StreetType %in% c("Highway", 
                            "Highway Interchange",
                            "Ramp")) -> highways

And Plotting

Visualizing A Single Entity

Using Normal Filtering

roads %>% 
  filter( FullName == "Three Chopt Road") %>%
  ggplot() + 
  geom_sf( aes(color=SpeedLimit), 
           lwd=2) + 
  coord_sf()

Intersection Plotting

Let’s grab a bit of zoning from The Fan

districts %>%
  filter( OBJECTID == 368 ) %>%
  st_buffer(dist = 1500) %>%
  st_bbox() -> fan_bbox

districts %>%
  st_crop( fan_bbox ) -> theFan 

plot( theFan["Name"])

Add Auxillary Data

zone_url <- "https://raw.githubusercontent.com/dyerlab/ENVS-Lectures/master/data/DistrictCodes.csv"

theFan %>%
  left_join( read_csv( zone_url ),
             by="Name") %>%
  mutate( Category = factor( Category) ) %>%
  select( OBJECTID, 
          Name, 
          Category, 
          everything() )  -> theFan

 

ggplot( theFan ) + geom_sf( aes( fill=Category)) + 
  scale_fill_brewer( type="qual", palette = "Set3")

Cropping Roads

roads %>% st_crop( st_bbox( theFan )) -> fanRoads
plot( st_geometry( fanRoads ))

Selecting a Single Zone

Add an attribute to the POLYGON’s indicating if the district is that big one in the middle of the plot.

theFan %>%
  mutate( Target = ifelse( OBJECTID == 368, 
                           TRUE, 
                           FALSE) ) -> theFan

Plot it

theFan %>%
  ggplot() + 
  geom_sf( aes(fill=Target) ) + 
  geom_sf_text( aes(label=OBJECTID), size=3 ) +
  coord_sf() + theme(legend.position = "none")

Isolate the District

target <- theFan[ theFan$OBJECTID == 368, ]
plot( st_geometry( target ) )

 

fanRoads %>%
  filter( st_intersects( fanRoads, target, 
                         sparse = FALSE ) == TRUE ) %>%
  as_data_frame() %>% select( `Street Name` = FullName ) %>%
  arrange( `Street Name`) %>% unique() 
# A tibble: 39 × 1
   `Street Name` 
   <chr>         
 1 Allison St    
 2 Birch St      
 3 Boyd St       
 4 Floyd Ave     
 5 Grove Ave     
 6 Hanover Ave   
 7 Kensington Ave
 8 Lombardy Pl   
 9 Madumbie Lane 
10 Monument Ave  
# … with 29 more rows

 

fanRoads %>%
  filter( st_intersects( fanRoads, target, 
                         sparse = FALSE ) == TRUE ) %>%
  ggplot() +
  geom_sf( data=target  ) +
  geom_sf( color="darkgreen" )

Questions

If you have any questions, please feel free to either post them as an “Issue” on your copy of this GitHub Repository, post to the Canvas discussion board for the class, or drop me an email.

Peter Sellers looking bored

Dyerlab